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The  nozzle  simulation  code  developed  for  AEDC  is  designed  for  the  analysis  of  hyper¬ 
sonic  converging-diverging  nozzles  including  the  effects  of  high  pressure,  finite-rate  chemical 
reactions  and  vibrational  relaxation,  and  turbulent  wall  boundary  layers.  In  addition,  a 
module  has  been  written  to  predict  the  onset  of  water  vapor  condensation  for  hydrocarbon 
combustion  product  test  gases.  In  the  following  section,  the  physical  models  are  discussed; 
additional  details  will  be  available  in  a  paper  under  preparation  for  publication. 


1.1  High-Pressure  Model 

The  nozzle  simulation  code  uses  the  excluded  volume  thermal  equation  of  state  (also 
known  as  the  Abel  or  Clausius  equation  of  state).  This  is  the  first  correction  to  the  perfect 
gas  equation  of  state,  and  it  accounts  for  the  volume  occupied  by  the  gas  particles.  The 
thermal  equation  of  state  is  written  as1 

pRT 

P  ~  1  -  Kp/M 

where  R  is  the  mixture  gas  constant,  ba  is  the  co-volume  of  the  mixture,  and  M  is  the 
mixture  average  molecular  weight.  For  a  reacting  gas,  b0/M  is  a  constant,  making  the 
model  depend  on  only  a  single  parameter.  b0  is  obtained  from  the  effective  molecular 
diameter,  a,  with 

b0  =  tv  a3  N0 

where  N0  is  the  Avogadro  number.  The  model  works  well  for  high  pressure  flows  that  are 
not  near  the  critical  state,  as  is  the  case  for  hypersonic  nozzle  flows. 

Note  that  the  excluded-volume  equation  of  state  must  be  used  when  computing  the 
nozzle  supply  conditions.  In  addition,  if  the  specific  enthalpy  of  the  gas  mixture  is  required, 
the  correct  expression  for  enthalpy  must  be  used: 


h  —  e  +  p/p  —  e  + 


1  -  b0p/ M 


The  code  includes  the  high-pressure  correction  in  the  inflow  boundary  condition. 

The  value  of  the  parameter  ba/M  is  obtained  from  fitting  to  high-accuracy  equa¬ 
tion  of  state  data  over  the  range  of  interest.2,3  For  pure  N2  the  appropriate  value  is 
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ba/M  =  0.001120  m3/kg  and  for  Air  ba/M  —  0.001024  m3/kg.  For  other  gases,  the  ef¬ 
fective  molecular  diameter  may  be  obtained  from  Table  1  or  other  sources.  The  mixture 
value  is  then  computed  by  number-weighting  (or  mole- weighting)  the  pure  gas  values: 


K 

M 


=  \*N,  ^ 

mix  O 


2^ 


<n_ 

Mo 


where  the  sum  is  over  all  species  in  the  gas  mixture,  and  js  is  the  mole  fraction  of  species 

s. 


Gas 

Symbol 

Molecular 
Diameter  (A) 

Argon 

Ar 

3.8 

Carbon  Dioxide 

co2 

2.8 

Carbon  Monoxide 

CO 

2.8 

Ethane 

c2H6 

4.4 

Helium 

He 

2.0 

Hydrogen 

h2 

2.4 

Methane 

ch4 

4.4 

n-Butane 

c4h10 

4.9 

Nitrogen 

n2 

3.0 

Oxygen 

02 

2.8 

Propane 

c3H8 

4.9 

Water 

h2o 

2.8 

Table  1.  Effective  molecular  diameters  for  relevant  gases;  from  www.domnickhunter.com. 


1.2  Internal  Energy 

In  hypersonic  nozzle  flows,  the  vibrational  energy  modes  of  the  gas  are  often  excited. 
The  nozzle  code  models  the  vibrational  excitation  using  a  simple  harmonic  oscillator,  which 
is  valid  for  typical  reservoir  temperatures.  The  gas  mixture  vibrational  energy  per  unit 
volume,  Ev.  is  given  by 


S 


Ps&vs 


where  the  vibrational  energy  for  each  species  is 


R 


Or 


Ms  2-;9rsexp (9vrs/Tv)  -  1 
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where  r  is  summed  over  each  vibrational  mode  of  species  s,  9vrs  is  the  characteristic 
temperature  of  vibration  for  that  mode,  and  grs  is  the  degeneracy  of  that  mode.  These 
constants  are  well  known  and  are  available  from  sources  such  as  the  NIST-JANAF  Thermo- 
chemical  Tables4  (previously  the  JANAF  Thermochemical  Tables).  In  the  above  expres¬ 
sion,  Tv  is  the  vibrational  temperature  that  characterizes  the  gas  mixture.  Thus,  we  have 
assumed  that  the  vibrational  modes  of  the  gas  relax  with  one  another  relatively  quickly, 
resulting  in  a  single  vibrational  temperature. 

The  vibrational  energy  modes  are  assumed  to  relax  toward  the  translational  modes 
using  the  Landau- Teller  model.  We  compute  a  number- weighted  relaxation  time  using  the 
individual  species  relaxation  times  from  Millikan  and  White.5  This  model  is  accurate  for 
most  species  pairs,  however  it  is  known  to  be  inadequate  for  carbon  dioxide  relaxation 
because  of  resonant  energy  transfer  between  modes.  Therefore,  we  use  the  results  of 
Camac6  to  improve  the  vibrational  energy  relaxation  model  when  carbon  dioxide  is  present. 

To  construct  the  total  energy  per  unit  volume  of  the  mixture,  E,  we  sum  the 
translational-rotational  energy,  the  vibrational  energy,  and  the  kinetic  energy 

E  —  pscvsT  +  Ev  +  —  p(u2  +  v2)  +  y]  psh° 

S  S 

where  cvs  is  the  translational-rotational  specific  heat  at  constant  volume  for  species  s.  For 
example,  diatomic  molecules  have  cV8  —  §-^-,  and  atoms  cvs  =  h°  is  the  heat  of 

formation  of  species  s. 

The  code  is  set  up  so  that  the  molecular  parameters  are  input  variables,  allowing  new 
gas  models  to  be  constructed  relatively  easily. 


1.3  Finite-Rate  Chemical  Reactions 


The  chemical  composition  of  the  gas  changes  as  it  flows  through  the  converging- 
diverging  nozzle.  In  many  cases,  the  expansion  is  so  sudden  that  an  chemical  equilibrium 
assumption  is  not  valid,  and  we  must  consider  finite-rate  recombination  of  the  test  gas.  To 
do  so,  we  solve  a  mass  conservation  equation  for  each  chemical  species,  including  a  mass 


production  or  destruction  term  due  to  chemical  reactions.  The  equation  is 


dps  _d_ 

dt  dxj 


(yPSUj  T  PsVsj) 


=  Ws 
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where  we  have  used  index  notation.  Here  vSj  is  the  diffusion  velocity  of  species  s  relative 
to  the  mass-averaged  velocity,  Uj.  ws  is  the  chemical  source  term  due  to  reactions;  this  is 
constructed  from  the  law  of  mass  action  and  an  appropriate  chemical  kinetics  model. 

For  typical  hydrocarbon  reaction  products,  we  consider  an  11-species  CHON  model 
consisting  of  C02,  H20,  CO,  N2,  02,  H2,  OH,  O,  H,  and  N.  The  kinetics  model  used  is 
based  on  rates  taken  from  the  Gas  Research  Institute  Mechanism  3.0  (GRI-Mech3.0).  The 
reactions  considered  are: 


H2+M^H  +  H  +  M 
H20  +  M  ^  OH  +  H  +  M 
OH  +  M^O  +  H  +  M 
02+M^0  +  0  +  M 
C02  +  M  ^  CO  +  O  +  M 
NO  +  M^N  +  O  +  M 
N2+M^N  +  N  +  M 
02  +  H  ^  OH  +  O 
H2  +  O  ^  OH  +  H 
OH  +  OH  ^  H20  +  H 
OH  +  H2  ^  H20  +  H 
OH  +  CO  ^  C02  +H 
NO  +  N  ^  N2  +  O 
02  +  N  ^  NO  +  O 
02  +  C0^0  +  C02 
N  +  OH  ^  NO  +  H 


(Rl) 

(R2) 

(R3) 

(R4) 

(R5) 

(R6) 

(R7) 

(R8) 

(R9) 

(RIO) 

(Rn) 

(R12) 

(R13) 

(R14) 

(R15) 

(R16) 


For  air  flows,  we  use  a  5-species  model  with  5  reactions  (R4,  R6,  R7,  R13  and  R14); 
and  for  a  pure  nitrogen  flow  we  just  use  the  nitrogen  dissociation/recombination  reaction 
(R7). 
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The  forward  reaction  rates  for  each  of  these  reactions  is  from  GRI-Mech3.0,  and  we 
compute  the  backward  rates  using  the  forward  rate  and  the  equilibrium  constant.  We  curve 
fit  the  equilibrium  constant  using  a  five-parameter  fit  over  the  temperature  range  of  interest 
by  minimizing  the  Gibbs  free  energy  based  on  the  Gordon  and  McBride  thermochemical 
data.7  The  curve  fits  provided  with  the  code  are  valid  in  the  range  200  <  T  <  3200 K; 
a  different  temperature  range  can  be  fit  with  minimal  effort  as  needed.  All  equilibrium 
constants  have  been  verified  by  comparison  with  the  GRIMech3.0  tabulations. 

1.4  Transport  Properties  and  Turbulence  Modeling 

The  molecular  viscosity,  thermal  conductivity,  and  mass  diffusivity  are  computed  using 

the  standard  Wilke  mixing  rule8,  with  pure  species  data  obtained  from  the  Blottner  curve 

fits.9  All  relevant  data  are  provided  in  the  input  files  for  the  nozzle  simulation  code.  In 

the  case  of  nitrogen  gas,  the  viscosity  fits  from  AEDC  are  used 

A  2.8676  x  10“8  T0-979  T  <  227  R; 

At  =  <  6.8971  x  10-7  T1-5/(179.1  +T)  227  <T<  795  R; 

l  1.9466  x  10“7  T0-659  T  >  795  R. 

where  T  is  in  Rankine  and  n  is  in  lbm/fts  (multiply  by  1.4882  to  obtain  SI  units  of  kg/m  s). 

The  translational-rotational  thermal  conductivity  is  obtained  through  the  use  of  an 

Eucken  relation,  which  is  essentially  a  variation  of  the  Prandtl  number  to  account  for  the 

mixture  of  different  components.  The  thermal  conductivity  of  pure  species  s  is 

/9  x  R 

ks  =  A4n^  +  °vs) 

The  mixture  value  is  then  obtained  using  the  Wilke  mixing  approach.  The  vibrational 
energy  conduction  is  written  as: 

Edevs 

rlvsl-ls  q  x 
S  3 

which  is  more  accurate  and  stable  than  basing  the  vibrational  heat  flux  on  the  vibrational 
temperature.  The  value  of  T)vs  is  taken  to  be  1.2  from  Ref.  10.  The  mass  diffusivity  is 
based  on  simple  Fickian  diffusion  with  a  constant  Lewis  or  Schmidt  number. 

It  is  critical  to  predict  the  correct  boundary  layer  thickness  in  hypersonic  nozzles  be¬ 
cause  it  determines  the  effective  area  ratio  of  the  nozzle,  and  therefore  the  facility  operating 
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conditions.  In  hypersonic  facilities,  the  nozzle  wall  boundary  layers  are  almost  always  tur¬ 
bulent  because  of  the  extremely  high  unit  Reynolds  number  in  the  throat  region.  However, 
conventional  algebraic  turbulence  models  ( e.g .  Baldwin-Lomax)  do  not  work  well  because 
of  the  highly  favorable  pressure  gradient  and  strong  temperature  gradients  in  the  boundary 
layer.  We  have  performed  extensive  comparisons  between  pitot  pressure  rake  measurements 
and  nozzle  wall  heat  transfer  rate  measurements  in  the  CUBRC  LENS  facilities.11  We  find 
that  the  one-equation  Spalart-Allmaras  turbulence  model12  with  the  Catris  and  Aupoix 
compressibility  correction13  gives  good  agreement  with  experiments  across  a  wide  range 
of  nozzle  operating  conditions.  This  is  the  turbulence  model  that  is  implemented  in  the 
APTU  nozzle  code. 

The  code  initializes  the  inflow  turbulent  viscosity  according  to  the  criteria  given  in 
Ref.  13.  There  are  no  other  adjustable  constants  in  the  present  implementation;  the 
turbulent  boundary  layer  is  allowed  to  develop  along  the  converging  section  of  the  nozzle 
according  to  the  local  flow  conditions.  Again,  because  hypersonic  nozzles  operate  at  high 
pressure  stagnation  conditions,  the  unit  Reynolds  numbers  are  so  high  that  the  boundary 
layer  is  fully  developed  before  the  flow  reaches  the  nozzle  throat. 

1.5  Water  Vapor  Condensation 

When  hydrocarbon  combustion  products  are  expanded  through  a  nozzle,  the  water 
vapor  may  condense  into  droplets  which  cause  a  sudden  entropy  rise  and  adversely  affect 
tunnel  operation.  We  have  implemented  a  condensation  onset  prediction  model  in  the  code 
based  on  the  work  of  Erickson  et  al.14 

The  main  criterion  for  condensation  onset  is  based  on  the  local  saturation  temperature 
of  the  gas  mixture.  Once  the  gas  becomes  saturated,  nucleation  will  begin.  Under  most 
hypersonic  nozzle  flow  conditions,  the  nucleation  rate  is  very  large  and  the  gas  condenses 
suddenly  when  it  becomes  saturated.14  Therefore,  if  we  monitor  the  saturation  temperature 
relative  to  the  gas  temperature,  we  can  determine  where  condensation  if  likely  to  occur.  In 
addition,  Ref.  14  gives  a  model  for  the  nucleation  rate,  J,  which  can  be  used  to  estimate 
the  rate  at  which  condensation  will  occur.  However,  in  practice  this  rate  is  typically  so 
large  that  condensation  is  very  rapid  once  the  gas  is  saturated.  (Particularly  because 
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the  gas  continues  to  expand  in  the  nozzle,  making  the  flow  more  and  more  susceptible  to 
nucleation  and  condensation.) 

We  have  implemented  the  Erickson  et  al.  model  for  a  combustion  gas  mixture.  The 
saturation  temperature,  Tsat  is  found  by  solving  the  equation 

6641.7 

£»h2o  =  exp (55.897 - — -  4.4864  lnTsat) 

J-  sat 


where  pn2o  is  the  local  partial  pressure  of  water  vapor  in  the  flow.  This  expression  comes 
from  a  curve-fit  to  the  vapor  pressure  of  water  at  a  given  temperature.  Therefore,  when 
Tsat  =  T,  the  gas  is  saturated  and  can  begin  to  nucleate  and  condense.  For  Tsat  >  T,  the 
gas  is  supersaturated  and  nucleation  will  be  very  rapid. 

The  nucleation  rate  can  be  found  by  first  computing  the  critical  droplet  radius,  r*, 
for  water  vapor 

2cru, 

v  — - 

Pi  -RT 1 2  O  111  (p>!  1 2  o  / pv  ) 

where  pe  is  the  density  of  liquid  water  pu20  is  the  partial  pressure  of  water  vapor,  and  pv 
is  the  vapor  pressure  of  water  over  a  flat  surface.  An  expression  for  pv  was  fit  in  Ref.  14  as 

6641.7 

pv  —  exp (55.897 - — - 4.4864  lnT). 


In  the  above  expression,  aw  is  the  surface  tension  of  water  given  by 


(7VJ  =  (82.27  +  75.612  TR  -  256.889  TR  +  95.928  TR)  x  1(T3 


where  TR  —  T/Tc  with  Tc  =  647.3  K  the  critical  temperature  of  water.  Droplets  larger 
than  r*  tend  to  grow,  while  droplets  less  than  r*  evaporate. 

Erickson  et  al.  derive  an  expression  for  the  nucleation  rate  of  water  in  a  carrier  gas. 
Their  expression  is 


J 


Qc  pLo  I  2aw 

1 +  Q  Pi  V^h.o6^1  3fcT  j 


where  qc  is  the  condensation  coefficient  (taken  to  be  unity  in  Ref.  14)  and 


2(7h2o  -  1)  L  /  L  _  U 
7h2o  +  1  -Rh2o^  Rr2oT  2 
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is  the  nonisothermal  correction  factor.  L  is  the  latent  heat  of  evaporation  at  T 

L  =  hn2o  ~  4.2  x  103  T  +  17.1175  x  106 

where  /ih2o  is  the  enthalpy  of  gaseous  water  and  the  units  of  L  are  J/kg.  Therefore,  we 
can  evaluate  the  rate  of  formation  of  water  droplets  as  a  function  of  the  local  state  of  the 
gas  in  the  nozzle  expansion.  A  sudden  increase  in  J  indicates  onset  of  condensation  and 
rapid  formation  of  water  droplets. 

Note  that  we  do  not  need  to  compute  any  of  these  quantities  unless  Tsat  >  T.  Also, 
we  do  not  couple  the  nucleation  and  droplet  growth  processes  to  the  flow,  and  therefore 
this  calculation  can  be  carried  out  as  a  post-processing  step. 

1.6  Boundary  Conditions 

The  nozzle  wall  boundary  conditions  are  straight-forward:  no-slip  for  velocity  and 
either  adiabatic  or  isothermal.  In  the  isothermal  case,  either  the  wall  temperature  can  be 
set  to  a  single  value,  or  a  list  of  wall  temperatures  at  each  grid  point  may  be  read  in. 

The  inflow  boundary  condition  is  more  difficult  to  set  up  because  of  vibrational  ex¬ 
citation  and  the  high  pressure  equation  of  state.  The  original  version  of  the  code  used 
a  characteristic-based  inflow  boundary  condition,  however  it  was  difficult  to  obtain  the 
correct  stagnation  conditions  from  the  inflow  plane  into  the  solution  domain. 

The  current  version  of  the  code  uses  a  different  approach.  The  total  mass  flux  across 
a  plane  of  the  interior  solution  is  computed,  and  the  inflow  velocity  is  then  set  so  that  the 
inflow  has  the  correct  mass  flux;  the  temperature  is  then  adjusted  to  maintain  the  specified 
stagnation  enthalpy  of  the  inflow  based  on  assuming  that  the  density  varies  only  weakly 
as  the  gas  accelerates.  This  approach  then  assures  that  the  inflow  conditions  are  correctly 
specified  and  that  once  the  nozzle  has  choked,  the  mass  flow  through  the  inlet  is  correct. 
Extensive  testing  has  validated  this  approach. 
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2.  Nozzle  Simulation  Tools 

The  nozzle  analysis  tools  comprise  a  computational  fluid  dynamics  solver,  codes  to 
aid  in  determining  reservoir  conditions,  a  grid  generator,  and  a  code  to  generate  plotting 
files.  In  this  section,  code  input/output  files  are  discussed,  specific  compilation  instructions 
for  various  computers  are  given,  and  the  steps  for  obtaining  a  solution  are  presented.  In 
the  next  section  several  test  cases  are  documented,  with  specific  instructions  to  set  up  the 
codes,  obtain  reservoir  conditions,  run  cases,  and  plot  results.  SI  units  are  used  throughout 
the  codes  except  where  noted.  The  nozzle  analysis  codes  are: 

•  aptu mozzle  .  f  CFD  code  for  hypersonic  nozzle  simulation  and  analysis.  The  code  is 
set  up  to  run  either  N2,  Air  (5-species  or  non-reacting),  and  carbon- hydrogen-oxygen- 
nitrogen  (CHON)  combustion  product  gas  mixtures. 

•  n2eq.f ,  aireq.f ,  choneq.f  Codes  to  compute  the  equilibrium  gas  composition  of 
nitrogen,  air,  or  CHON  given  a  pressure,  temperature  and  elemental  mixture.  It  uses 
the  same  curve-fits  for  the  equilibrium  constants  as  the  nozzle  code,  so  that  the  inflow 
conditions  are  compatible.  The  codes  use  the  excluded  volume  model  for  high  pressure 
flows  used  in  the  CFD  code. 

•  read.f  Code  to  read  solution  files  and  write  out  a  Tecplot-compatible  data  file.  For 
the  hydrocarbon  combustion  product  simulations,  the  saturation  temperature  and 
nucleation  rate  are  written  out. 

•  nozgrid.f  This  is  a  simple  elliptical  grid  generator  that  takes  as  input  a  list  of  x 
and  y  nozzle  coordinates  and  generates  a  nozzle  grid.  It  adds  a  constant-area  section 
upstream  of  the  converging  section,  which  works  better  with  the  inflow  boundary 
condition. 
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2.1  Grid  Generation 

To  run  the  CFD  code,  you  first  need  to  generate  a  grid  using  nozgrid .  f  or  some  other 
grid  generator,  nozgrid. f  expects  to  read  a  datafile  of  x  and  y  data  (in  inches)  of  the 
form  given  in  the  file  coords.dat.*  Note  that  the  last  line  of  coords.dat  corresponds  to 
-999  -999  -999  to  signify  the  end  of  the  data.  It  is  up  to  the  user  to  provide  a  reasonable 
axial  variation  of  the  x-y  pairs.  In  general,  the  nozzle  coordinates  should  be  clustered  near 
the  throat  and  in  the  regions  of  rapid  slope  change.  A  good  nozzle  grid  will  typically  have 
about  500  x-y  pairs;  the  number  of  grid  points  typically  increases  with  the  operating  Mach 
number  of  the  nozzle. 

Within  nozgrid. f,  the  wall  normal  spacing  can  be  adjusted  by  the  user  by  changing 
the  value  of  rk;  larger  values  increase  the  grid  stretching  (decrease  the  wall-normal  spac¬ 
ing).  The  elliptical  grid  generator  is  very  sensitive  to  the  value  of  this  parameter,  so  only 
small  changes  should  be  made  to  rk. 

Running  nozgrid. f  results  in  the  hie  nozzle. grid  of  the  correct  form  for  use  with 

the  CFD  solver.  Note  that  the  CFD  method  uses  dummy  (or  ghost)  cells  to  enforce  the 

boundary  conditions,  nozgrid. f  generates  a  grid  that  includes  the  dummy  cells  and  writes 

the  grid  according  to  the  following  format: 

open (21  ,file=;’ nozzle  .  grid’ ) 
write (21,*)  il 
write (21,*)  jl 

write (21,*)  ((xg(j , i) , i=l,il+l) , j=l, jl+1) , 

( (yg ( j ,i> ,i=l,il+l) , j=l, jl+1) 

close (21) 

Where  il  and  jl  are  the  number  of  finite  volume  cells  (including  the  dummy  cells)  in 
the  axial  and  normal  directions,  respectively.  If  a  different  grid  generator  is  used,  the  grid 
must  be  written  in  this  form  including  the  dummy  cells.  To  form  the  dummy  cells,  simply 
mirror  the  grid  points  about  the  inflow/outflow  planes,  nozzle  centerline,  and  nozzle  wall; 
the  values  of  il  and  jl  then  to  need  to  be  increased  appropriately. 

*  The  format  is  a  list  of  three  real  numbers  per  line:  x,  y.  and  another  quantity  that  is 
not  used.  The  units  are  inches  for  x  and  y\  a  straight-forward  edit  of  nozgrid .  f  will  allow 
conversion  from  other  units. 
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Please  note  that  the  grid  generated  by  nozgrid.f  must  be  checked  by  the  user  to 
verify  that  the  grid  generation  has  been  performed  correctly.  Given  the  wide  range  of 
hypersonic  nozzle  geometries,  it  is  very  difficult  to  develop  a  grid  generation  code  that  will 
work  well  for  all  cases. 

2.2  Stagnation  (Reservoir)  Conditions 

Codes  are  provided  in  the  distribution  to  allow  the  user  to  compute  the  appropriate 
stagnation  conditions  for  high-pressure  mixtures  of  N2,  Air,  and  CHON  mixtures.  These 
codes  (n2eq.f,  aireq.f,  and  choneq.f)  determine  the  equilibrium  composition  of  the 
gas  by  solving  the  polynomial  system  of  equations  resulting  from  the  definition  of  the 
equilibrium  constants.  The  user  must  input  the  elemental  composition  of  the  gas  for  the 
CHON  mixture;  this  is  obtained  from  the  known  fuel  composition  and  fuel/air  mixture  of 
the  combustion  products.  The  solver  then  determines  the  equilibrium  composition  at  the 
input  stagnation  pressure  and  temperature.  The  output  from  these  codes  is  then  used  in 
the  input  file  for  the  nozzle  analysis  code.* 

The  CHON  equilibrium  solver  uses  the  LAPACK  linear  algebra  routines,  and  therefore 
it  must  be  linked  to  them  at  compile.  For  Portland  Group  Fortran90  use: 
pgf90  -r8  -fast  -o  choneq.x  choneq.f  -llapack  -lblas 

2.3  Code  Compile 

To  compile  the  CFD  code,  you  must  use  MPI  and  Fortran90,  and  you  want  to  use 
whatever  optimization  is  available.  Double  precision  must  be  specified  during  the  compile. 
The  nozzle  analysis  code  has  been  tested  on  a  variety  of  machines  and  compilers: 

For  Portland  Group  (PGI)  Fortran  on  a  Linux  (Pentium  or  Opteron)  cluster: 

mpif90  -r8  -fast  -o  aptumozzle  .x  aptu mozzle .  f 
For  Intel  Fortran  on  a  Linux  cluster: 

*  Note  that  because  the  equilibrium  codes  solve  for  the  roots  of  the  polynomials  in 
mixture  mole  fraction,  sometimes  spurious  roots  of  these  equations  may  be  found.  The 
codes  have  been  tested  at  a  range  of  conditions,  but  the  user  must  be  aware  of  this  issue; 
improved  initial  estimates  of  the  species  mole  fractions  will  overcome  this  problem. 
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mpif90  -r8  -axW  -o  aptumozzle  .x  aptumozzle .  f 
The  -axW  flag  can  be  adjusted  depending  on  the  processor  on  the  machine. 

For  Intel  Fortran  on  an  SGI  Altix: 

ifort  -r8  -03  -o  aptumozzle.  x  aptumozzle.  f  -lmpi 
More  aggressive  optimization  should  work  too. 

2.4  Input  File 

Next,  you  need  to  specify  the  run  conditions;  this  is  done  in  the  file  nozzle,  inp. 
First,  let  us  focus  on  setting  the  reservoir  conditions  in  the  input  file.  Use  the  example 
input  files  from  the  nitrogen,  air,  and  CHON  mixtures  examples  as  a  starting  point.  This 
will  reduce  the  difficulty  of  providing  the  necessary  species  data  for  each  mixture.  The 
variables  in  nozzle .  inp  that  need  to  be  changed  are: 

density:  stagnation  density  in  kg/m3 

Tin,  Tvin:  stagnation  temperature,  vibrational  temperature  in  K 

cs :  inflow  mass  fraction  for  each  species  in  the  order  given  in  the  relevant  input  file  (N2, 
air,  or  CHON). 

iwall :  for  an  adiabatic  nozzle  wall  (seldom  relevant),  set  iwall  =  1.  For  a  single  uniform 
temperature  isothermal  nozzle  wall  set  iwall  =  2,  and  set  twall  to  the  appropri¬ 
ate  wall  temperature  in  Kelvins.  And  for  a  variable  temperature  isothermal  wall, 
set  iwall  =  3.  In  the  latter  case,  the  wall  temperature  values  (in  Kelvins)  must 
be  specified  in  an  external  file  (twall.dat  in  the  example  cases).  These  values 
correspond  to  the  wall  temperature  at  the  center  of  each  of  the  il  axial-direction 
finite  volumes.  The  name  of  this  file  is  specified  in  the  6th  line  of  the  input  file, 
iturb:  switch  to  turn  on/off  the  Spalart-Allmaras  turbulence  model.  Set  iturb  =  1  for 
a  turbulent  boundary  layer. 

2.5  Code  Execution 

Now  it  is  (finally)  time  to  run  the  nozzle  simulation  code.  The  CFL  number  sequence 
and  number  of  time  steps  are  the  main  variables  that  must  be  adjusted  during  a  run.  These 
should  be  set  so  as  to  obtain  a  converged  (steady-state)  solution  in  the  minimum  number 
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of  time  steps.  In  general,  the  solution  progresses  by  first  reaching  a  steady-state  in  the 
converging  section  by  establishing  the  choked  condition  at  the  throat.  Then  the  solution 
propagates  down  the  supersonic  portion  of  the  nozzle  until  the  solution  converges  and 
the  initial  condition  is  swept  out  of  the  computational  domain.  This  process  mimics  the 
physical  nozzle  starting  process  in  a  hypersonic  facility.  It  is  fundamentally  impossible  for 
the  nozzle  to  be  converged  until  a  physically  meaningful  time  has  elapsed  in  the  simulation. 
Therefore  it  is  advantageous  to  run  at  the  maximum  stable  time  step  to  reduce  the  number 
of  time  steps  to  obtain  a  fully  converged  (started)  nozzle  flow. 

However,  because  the  problem  is  highly  non-linear,  it  is  not  possible  to  simply  start 
the  simulation  at  an  arbitrarily  large  time  step.  Rather,  the  time  step  must  be  ramped 
up  slowly  so  that  the  solution  remains  physically  meaningful  (non-negative  density  and 
temperature)  during  the  initial  transient.  Unfortunately,  because  each  problem  is  different, 
the  optimal  ramping  schedule  varies  with  the  reservoir  conditions,  the  grid  size,  and  most 
importantly,  the  degree  of  grid  stretching.  The  examples  give  typical  time  step  ramping 
functions  for  several  cases. 

You  shouldn’t  need  to  change  most  of  the  switches  in  the  input .  inp  file;  you  will 
need  to  change  the  following  variables: 

istop:  number  of  time  steps  to  be  run  before  the  code  stops  and  writes  the  solution. 
The  number  of  time  steps  required  depends  on  the  problem.  As  mentioned  above, 
factors  that  affect  the  number  of  iterations  required  to  fully  start  the  nozzle  flow 
include  the  grid  size,  the  physical  length  of  the  nozzle,  the  area  ratio  of  the  nozzle, 
and  the  rate  at  which  the  user  ramps  up  the  CFL  number  (while  maintaining  a 
physically  meaningful  solution).  For  the  cases  I  have  run,  the  calculation  converges 
in  about  500  time  steps  for  a  short,  low  Mach  number  nozzle  to  3000  for  a  long  high 
Mach  number  nozzle. 

iconr:  switch  to  start  from  previous  solution  (iconr  =  1)  or  from  a  guessed  initial  solu¬ 
tion  (iconr  ^  1).  Generally  set  iconr  =  -1  and  run;  if  the  solution  needs  to  be 
run  farther,  set  iconr  =  1  and  run  some  more. 

Additional  information  on  the  run-time  parameters  is  given  in  nozzle,  inp. 
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The  ramping  function  for  the  CFL  number  sequence  is  given  following  the  species 
information  in  input .  inp.  This  list  tells  how  the  code  ramps  the  CFL  number  from  the 
initial  value  given  as  the  variable  cf  1.  Every  20  time  steps  the  code  takes  the  next  value 
of  CFL  in  the  sequence  up  to  the  point  where  it  detects  a  negative  value  of  CFL;  then 
it  continues  with  the  last  non-negative  value  until  the  run  is  complete.  Thus  in  a  typical 
case,  the  calculation  starts  at  0.1  and  after  20  iterations,  goes  to  1.0,  5.0,  up  to  1000.0 
or  usually  more.  The  exact  sequence  depends  on  the  problem  and  may  differ  from  what  I 
have  given  you.  If  you  can’t  run  at  cf  1  >  100  after  several  hundred  iterations,  there  is 
probably  something  wrong. 

From  experience,  the  maximum  CFL  number  is  usually  limited  by  a  physical  time  step 
size  -  typically  a  CFL  that  corresponds  to  about  a  time  step  of  several  microseconds.  Thus, 
for  highly  stretched  grids  with  very  small  wall-normal  spacing,  the  maximum  CFL  number 
can  be  large  (as  much  as  106) ,  while  for  lower  Reynolds  number  flows,  the  maximum  CFL 
may  reach  only  several  thousand. 

As  the  calculation  is  running,  several  things  should  be  monitored:  The  solution  resid¬ 
ual  (the  second  column  of  numbers  written  to  the  screen  and  to  converge. m)  should 
generally  trend  downwards,  though  it  may  go  through  periods  where  it  increases.  In  most 
cases,  the  residual  will  drop  by  8  or  more  orders  of  magnitude  as  the  solution  converges.  If 
the  residual  begins  to  increase  significantly  over  a  large  number  of  time  steps,  the  solution 
is  probably  diverging  and  the  CFL  is  too  large  or  has  been  increased  too  rapidly;  restart 
the  calculation  with  a  less  aggressive  CFL  ramping  function.  In  addition  to  the  residual, 
the  time  step  size  should  be  monitored  as  well  (the  third  column  of  numbers  written  to 
the  screen).  Again,  unless  the  calculation  has  progressed  though  a  physically  meaningful 
time,  it  cannot  be  converged  to  a  steady-state.  A  simple  estimate  of  the  nozzle  starting 
time  will  give  guidance  as  to  the  number  of  time  steps  required  to  computationally  start 
the  nozzle. 

In  some  cases,  the  residual  will  not  decrease  by  a  huge  (8+  order  of  magnitude)  factor 
due  to  oscillations  induced  by  the  turbulence  model  source  terms.  Experience  shows  that 
these  oscillations  do  not  affect  the  solution  quality,  and  therefore  the  residual  does  not 
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give  a  full  picture  of  whether  the  solution  is  converged.  In  practice,  the  solution  at  one 
or  two  points  will  oscillate  between  two  states,  resulting  in  a  relatively  large  residual  but 
no  adverse  effect  to  the  overall  solution.  That  is  why  it  is  important  to  also  monitor  the 
total  flow  time  to  determine  whether  the  nozzle  has  started  based  on  physical  arguments. 
The  new  code  has  been  significantly  improved  in  this  regard;  in  all  test  cases  performed 
to  date,  this  problem  has  not  occurred. 

Finally,  it  is  now  time  to  run  the  calculation!  Depending  on  how  your  system  does 
MPI,  you  should  be  able  to  do  something  like  (for  8  processors): 
mpirun  -np  8  aptujnozzle .x 

This  will  result  in  an  output  file  for  each  processor  of  the  name  nozzle  .  flow .  *,  which  can 
be  read  by  read.f. 

2.6  Post-processing  and  Plotting 

The  code  read.f  reads  the  output  files  from  the  CFD  code  and  writes  out  data  files  in 
Tecplot  compatible  format  for  plotting  purposes.  The  code  reads  nozzle .  inp  to  determine 
the  gas  model  (iset)  and  writes  data  appropriate  for  the  gas  model.  Again,  the  data  are 
written  in  SI  units. 

To  post-process  the  data,  compile  the  code  using  a  Fortran90  compiler.  For  Portland 
Group  Fortran90: 

pgf90  -r8  -fast  -o  read.x  read.f 

Then  simply  run  the  executable.  The  user  will  need  to  input  the  number  of  processors 
used  iu  the  simulation;  all  other  information  is  obtained  from  the  data  hies. 

For  the  CHON  mixture,  read.f  computes  the  saturation  temperature  and  nucleation 
rate  as  described  above.  If  additional  variables  need  to  be  plotted,  the  code  may  be 
modified  to  compute  them. 
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3.  Examples 

Included  with  the  distribution  are  4  example  problems:  N2,  Air,  and  two  CHON 
hypersonic  nozzle  flows.  This  section  describes  how  the  input  data  are  generated  for  each 
case,  and  details  of  each  run  are  provided. 

The  files  are  located  in  separate  directories:  N2example,  AIRexample,  HTTexample, 
and  CHONexample.  In  each  directory,  there  is  a  copy  of  the  grid  files,  input  files,  and 
solution  files.  The  same  version  of  aptu jnozzle . f  and  read.f  was  used  to  compute  each 
flow  field. 

3.1  Nitrogen  Example:  AEDC  Tunnel  9  Nozzle 

In  the  code  distribution,  this  example  problem  is  provided  in  the  directory  N2example. 
In  this  case,  a  301  x  121  (including  dummy  cells)  grid  was  provided  by  Greg  Molvik,  along 
with  the  wall  temperature  at  various  locations  along  the  nozzle.  A  small  code  (twall.f) 
was  written  to  linearly  interpolate  the  wall  temperature  for  each  of  the  wall  grid  cells. 

For  this  case,  the  stagnation  conditions  are:  Ta  =  3027  R  =  1681.5  K  and  pQ  = 
21,540psi  =  148.46  MPa.  Using  ba/M  =  0.001120  m3/kg,  the  density  is  computed  to  be 
223.145  kg/m3.  These  are  the  values  used  in  the  input  file  (nozzle .  inp). 

Because  the  nozzle  is  very  long  (>  12  m)  with  a  small  throat,  the  physical  nozzle  start 
time  is  expected  to  be  significant.  The  nominal  nozzle  exit  Mach  number  is  14,  which  for 
the  stagnation  temperature,  gives  a  test  section  axial  velocity  of  1850  m/s.  If  we  assume 
that  the  average  speed  of  the  gas  in  the  nozzle  is  half  of  this  value,  the  nozzle  starting  time 
is  approximately  12/925  =  0.013  sec.  Therefore,  our  calculation  must  be  run  for  at  least 
this  physical  time  before  we  can  expect  the  nozzle  to  be  started.* 

In  the  input  file,  there  is  a  fairly  aggressive  ramping  function  provided,  which  gives  a 
final  CFL  number  of  106,  resulting  in  a  time  step  of  approximately  10  /isec  =  10_5sec. 
Therefore,  once  we  have  reached  this  large  CFL,  we  expect  the  solution  to  require 
0.013/10-5  =  1300  time  steps.  With  a  factor  of  two  safety  factor,  we  run  the  solution 

*  The  time  integration  is  not  perfectly  time  accurate,  but  the  time  step  can  be  used  to 
estimate  the  elapsed  time  of  the  simulation.  In  practice,  it  has  been  found  that  it  does 
provide  useful  guidance  to  assess  the  progress  of  the  solution. 
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for  2500  time  steps. 

For  the  example  provided,  we  ran  on  4  dual-processor  Pentium  Xeon  2.4  GHz  machines 
connected  with  Myrinet.  The  compiler  was  Portland  Group  Version  4.1.  The  total  run 
time  was  550  seconds  for  2500  time  steps. 

The  output  files  are  provided  in  the  distribution,  including  the  output  from  read,  f .  A 
plot  of  the  L2  norm  of  the  solution  residual  is  provided,  which  is  from  the  file  converge  .m. 
Note  the  14  order  of  magnitude  decrease  in  the  residual  to  machine  zero  at  about  2000 
iterations.  In  this  case,  the  estimated  number  of  time  steps  was  reasonable  though  a  little 
optimistic. 


Figure  1.  Solution  residual  for  N2  test  case. 
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3.2  Air  Example:  Generic  5-Species  Air  Nozzle  Flow 

For  this  calculation,  we  use  the  same  Tunnel  9  nozzle  grid  that  was  used  for  the 
previous  example.  We  choose  typical  stagnation  conditions  of  pa  =  5496  psi  and  total 
enthalpy  hQ  —  2.9492  x  10'ft2/sec2.  These  conditions  are  inserted  into  aireq.f,  which 
is  compiled  and  run.  A  guess  for  the  stagnation  temperature  is  required  -  in  this  case 
something  like  3000  K  works.  The  code  then  writes  the  computed  stagnation  conditions 
to  the  screen,  which  are  then  copied  into  nozzle .  inp. 

For  this  case,  we  ran  a  constant  temperature  cold  wall  by  setting  iwall  =  2,  and 
twall  =  300.0.  In  addition,  iset  must  be  changed  to  2  to  use  the  5-species  air  model. 
The  CFD  code  is  then  run;  in  this  case,  we  ran  on  8  dual-processor  Xeons  with  a  total 
execution  time  of  547  sec.  The  increase  in  cost  is  due  to  the  5-species  air  model,  which 
involves  the  solution  of  more  equations  and  a  more  complex  chemical  kinetics  model. 

For  this  case  due  to  the  differences  in  the  stagnation  conditions  and  wall  temperature 
distribution,  the  time  step  that  corresponds  to  a  CFL  of  106  is  22  psec.  Therefore,  for 
this  run  the  maximum  CFL  was  reduced  by  a  factor  of  two.  In  this  case,  the  solutions 
does  not  converge  as  well  as  the  N2  example,  with  only  a  total  residual  reduction  of  8 
orders  of  magnitude.  This  does  not  necessarily  mean  that  the  solution  has  not  reached 
convergence;  in  this  case  the  chemistry  model  or  the  turbulence  model  is  causing  small 
levels  of  switching  between  solutions.  This  can  be  verified  by  restarting  the  simulation  and 
running  for  a  meaningful  number  of  time  steps  (say  200)  and  then  comparing  results. 
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3.3  CHON  Example:  NASA  Langley  8'HTT  Nozzle 

To  illustrate  and  validate  the  water  condensation  model,  we  simulated  a  case  from  the 
work  of  Erickson  et  al.14  for  the  NASA  Langley  S'  High  Temperature  Tunnel  nozzle.  Their 
approach  is  quasi-one-dimensional,  so  we  cannot  reproduce  their  results  exactly.  However, 
we  generated  a  nozzle  grid  from  the  information  provided  in  their  report.  Table  2  gives 
x-y  coordinates  of  the  diverging  section  of  the  nozzle  at  8  locations;  we  fitted  spline  curves 
through  these  points,  and  then  merged  the  converging  section  of  the  Tunnel  9  nozzle  to 
that  geometry.  Then  ran  nozgrid.f  to  produce  a  373  x  160  grid. 

The  conditions  for  the  run  are  provided  in  Ref.  14;  we  chose  to  consider  Case  1,  which 
has  Ta  =  1900  K  and  pQ  —  5.0  MPa.  This  case  has  methane-air  equivalence  ratio  of  0.798. 
To  obtain  the  relevant  mixture,  we  considered  the  stoichiometric  reaction: 

CH4  +  202  ->  2H20  +  C02 

For  an  equivalence  ratio  of  0.798  in  air,  we  then  have: 

0  79  0  79 

0.798CH4  +  202  +  — 2N2  -»■  0.798  x  2H20  +  0.798CO2  +  0.202  x  202  +  — 2N2 

U  •  ^  JL  U .  Zj  JL 

which  gives  a  mole  fraction  of  the  mixture  of:  7n2  =  0.728923,  7h2o  =  0.154624,  7co2  — 
0.077312,  and  7o2  =  0.039140.  This  mixture  is  put  into  choneq.f,  which  is  then  compiled 
and  run.  The  equilibrium  mixture  at  the  1900  K,  5  MPa  stagnation  conditions  is  then 
copied  into  input .  inp  and  the  variable  iset  =  3.  For  this  case,  we  do  not  know  the  wall 
temperature  distribution,  so  for  simplicity,  we  use  a  constant  temperature  isothermal  wall 
set  to  800  K. 

This  case  was  run  for  1000  iterations  at  a  maximum  CFL  of  4000.  The  relatively 
low  CFL  number  is  a  result  of  using  a  weakly  stretched  grid  due  to  the  lower  operating 
Mach  number  of  the  facility.  However,  this  CFL  results  in  a  75.5  //sec  time  step,  which 
rapidly  starts  the  nozzle.  In  this  case,  we  were  only  able  to  achieve  6  orders  of  magnitude 
convergence  because  of  grid  quality  issues  at  the  point-  where  the  constant  area  section 
joins  the  converging  section.  The  run  took  1048  seconds  on  16  processors. 

Figure  2  plots  contours  of  the  saturation  temperature  relative  to  the  local  temperature, 
Tsat  —  T;  the  gas  starts  to  condense  in  regions  where  this  quantity  is  positive.  The  figure 


20 


AEDC-TR-05-2 


also  plots  the  nucleation  rate,  J,  in  particles/m3  sec.  Clearly,  J  is  very  high  in  this  region 
and  the  gas  will  condense  very  rapidly.  These  results  show  Tsa t  >  T  at  the  same  location 
as  Ref.  14. 
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Figure  2.  Saturation  temperature  relative  to  the  local  temperature  (top)  and  nucleation 
rate  (bottom). 
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3.4  CHON  Example:  Mach  6.5  APTU  Nozzle 

The  final  example  is  for  the  Mach  6.5  APTU  nozzle  at  conditions  provided  by  Greg 
Molvik.  These  correspond  to  the  operating  condition  of  Ta  =  4800  R  and  pa  —  2800  psi. 
The  mole  fractions  are  obtained  by  running  choneq.f  starting  with  the  gas  composition 
provided.  Based  on  the  effective  molecular  diameters  (Table  1)  and  the  gas  composi¬ 
tion,  the  value  of  ba/M  is  computed  to  be  0.00076552  m3/kg.  This  gives  a  density  of 
25.533  kg/m3. 

The  nozzle  contour  data  file,  coords.dat,  has  a  poor  axial  distribution  of  x-y  coor¬ 
dinate  pairs.  Therefore,  a  short  code  was  written  to  read  in  the  data  file,  and  linearly 
interpolate  between  the  points  to  yield  a  better  distribution  of  grid  points.  This  code, 
nozcont .  f  is  provided  in  the  distribution.  This  code  is  first  run  to  generate  a  new  coor¬ 
dinates  file  (nozcont.dat),  and  then  the  grid  generation  file  nozgrid.f  is  run  to  generate 
the  grid.  For  the  parameter  settings  in  these  codes,  this  results  in  a  410  x  160  grid. 

The  CFD  code  was  run  for  1000  time  steps  ramping  to  a  maximum  CFL  of  2  x  105, 
corresponding  to  a  time  step  of  18  /isec.  The  residual  was  reduced  by  7  orders  of  magnitude, 
and  the  total  elapsed  physical  time  was  approximately  13  msec,  which  is  sufficient  to  fully 
converge  the  solution  for  this  relatively  short,  low  Mach  number  nozzle.  The  maximum 
residual  is  located  at  the  corner  where  the  constant  area  duct  suddenly  transitions  to  the 
converging  part  of  the  nozzle.  The  grid  in  this  region  is  highly  stretched,  which  likely  leads 
to  the  high  residual.  The  calculation  required  1136  sec  on  8  dual-processor  Xeons.  Figure 
3  plots  the  computed  Mach  number  contours  for  this  case. 

Figure  4  plots  the  solution  residual  for  two  grids.  The  poor  result  was  obtained  using 
the  original  coords .  dat  file,  and  the  significantly  improved  convergence  plot  was  obtained 
using  the  smoothly  varying  coordinates.  Clearly,  using  a  high  quality  grid  is  important 
in  obtaining  good  convergence  properties  and  a  rapid  solution.  Note  that  the  maximum 
allowable  CFL  on  the  poor  grid  was  just  1000. 
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Figure  3.  Mach  number  contours  for  the  Mach  6.5  APTU  nozzle  example. 


Figure  4.  Convergence  history  for  APTU  problem  on  two  different  grids. 
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4.  Code  Structure 

The  schematic  given  in  Fig.  5  illustrates  the  structure  of  the  nozzle  simulation  code. 
First,  several  subroutines  are  called  to  initialize  the  MPI  variables  (subroutine  startup), 
the  input  files  are  read  in  (subroutine  readin),  the  geometry  is  read  in  from  the  grid 
file,  and  the  mesh  metrics  are  computed  (subroutine  geometry),  and  finally  the  flow 
field  is  initialized  (subroutine  init).  When  the  code  is  starting  from  no  initial  solution, 
subroutine  nozzle_init  is  called  to  initialize  the  nozzle  flow  using  an  approximate  quasi- 
one- dimensional  solution  for  the  nozzle  geometry;  otherwise  the  solution  is  read  in  using 
subroutine  rdwrt. 

Then  the  solution  begins  by  running  a  series  of  istop  time  steps.  This  involves 
computing  the  inviscid  fluxes  across  each  finite  volume  surface;  this  is  done  in  subroutine 
fluxes.  This  subroutine  then  calls  the  viscous  flux  calculation,  the  chemical  source  term 
calculation,  and  then  the  Spalart-Allmaras  source  term  calculation.  The  inflow  boundary 
conditions  are  then  computed  in  the  subroutine  lbcuion. 

At  this  point,  all  required  information  for  a  time  step  has  been  computed.  The  time 
step  is  either  done  using  an  explicit  time  step  (iextst  =  1)  or  using  the  implicit  data- 
parallel  line-relaxation  method15  (iextst  ^  1).  In  the  latter  case,  the  subroutine  dplr  is 
called  to  assemble  the  Jacobian  matrices,  factorize  the  system  of  equations  (subroutine 
factor),  and  then  backward-substitute  into  the  block  tri-diagonal  system  (subroutine 
solve).  Details  of  the  linear  algebra  performed  by  these  routines  is  given  in  Ref.  15. 

Once  the  istop  time  steps  are  complete,  the  solution  is  written  out  using  subroutine 
rdwrt,  and  the  simulation  ends. 
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5.  Variable  Definitions 

The  following  is  a  list  of  the  major  variables  used  in  the  nozzle  simulation  code.  In  all 
cases,  the  variables  are  dimensional  and  use  SI  units;  these  are:  meters,  seconds,  Kelvins, 
Pascals,  Joules,  kilograms,  and  kilomoles. 


Grid  Variables 

X 

axial-direction  distance 

[m] 

y 

radial-direction  distance 

[m] 

si 

i-direction  surface  area  of  finite  volume 

H 

sj 

/-direction  surface  area  of  finite  volume 

H 

volin 

inverse  volume  of  finite  volume 

[m~2] 

Flow  Variables 

r 

density 

[kg/m3] 

rhos 

species  density 

[kg/m3] 

u 

axial-direction  velocity 

[m/s] 

V 

radial-direction  velocity 

[m/s] 

c 

speed  of  sound 

[m/s] 

t 

translational- rotational  temperature 

[K] 

tv 

vibrational  temperature 

[K] 

P 

pressure 

[Pa] 

re 

total  energy  per  unit  volume 

[J  /  m3] 

rev 

vibrational  energy  per  unit  volume 

[J  /  m3] 

evs 

species  vibrational  energy  per  unit  mass 

[J/kg] 

rmut 

turbulent  eddy  viscosity 

[kg/m  s] 

cv 

translational-rotational  specific  heat 

[J/kg  K] 

gcon 

mixture  gas  constant 

[J/kg  K] 

res 

L2  solution  norm  of  density 

[] 

twall 

wall  temperature 

[K\ 

Thermodynamic/Rate  Constants 

hf  orm 

heat  of  formation 

[J/kg] 
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thetv 

characteristic  temperature  of  vibration 

[K] 

vdg 

vibrational  degeneracy 

[] 

smw 

species  molecular  weight 

[kg/kmole] 

CVS 

species  translational-rotational  specific  heat 

[kg/kmole] 

rmu 

mixture  viscosity 

[kg/ms] 

rkap 

mixture  thermal  conductivity 

[J/ms  K] 

gsp 

species  gas  constant 

[J/kgK] 

cf 

forward  reaction  rate  coefficient 

[m3/kmoles 

aeta 

temperature  exponent  for  forward  reaction  rate 

[] 

athet 

exponential  constant  for  forward  reaction  rate 

[K\ 
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